
/* copyright (c) 1996, 1997, 1998, 1999 William R. Pearson and the
   U. of Virginia */

/* $Name: fasta34t11d4 $ - $Id: compacc.c,v 1.20 2002/01/03 18:31:29 wrp Exp $ */

/* Concurrent read version */

#include <stdio.h>
#include <stdlib.h>
#if defined(UNIX) || defined(MSDOS)
#include <sys/types.h>
#endif

#include <limits.h>

#include <string.h>
#include <time.h>
#include <math.h>

#include "defs.h"
#include "param.h"
#include "structs.h"

#ifndef PCOMPLIB
#include "mw.h"
#else
#include "p_mw.h"
#endif

/*
#define XTERNAL
#include "uascii.h"
#include "upam.h"
#undef XTERNAL
*/

#ifdef PCOMPLIB
#include "msg.h"
extern int nnodes;
#ifdef PVM_SRC
extern int pinums[],hosttid;
#endif
#ifdef MPI_SRC
#include "mpi.h"
#endif
#endif

extern time_t tdone, tstart;		/* Timing */
extern void abort ();
extern void ptime ();
/* extern void summ_stats(char *s_str); */

void selectbest(bptr,k,n)	/* k is rank in array */
     struct beststr **bptr;
     int k,n;
{
  int v, i, j, l, r;
  struct beststr *tmptr;

  l=0; r=n-1;

  while ( r > l ) {
    v = bptr[r]->score[0];
    i = l-1;
    j = r;
    do {
      while (bptr[++i]->score[0] > v) ;
      while (bptr[--j]->score[0] < v) ;
      tmptr = bptr[i]; bptr[i]=bptr[j]; bptr[j]=tmptr;
    } while (j > i);
    bptr[j]=bptr[i]; bptr[i]=bptr[r]; bptr[r]=tmptr;
    if (i>=k) r = i-1;
    if (i<=k) l = i+1;
  }
}

void selectbestz(bptr,k,n)	/* k is rank in array */
     struct beststr **bptr;
     int k,n;
{
  int i, j, l, r;
  struct beststr *tmptr;
  float v;

  l=0; r=n-1;

  while ( r > l ) {
    v = bptr[r]->zscore;
    i = l-1;
    j = r;
    do {
      while (bptr[++i]->zscore > v) ;
      while (bptr[--j]->zscore < v) ;
      tmptr = bptr[i]; bptr[i]=bptr[j]; bptr[j]=tmptr;
    } while (j > i);
    bptr[j]=bptr[i]; bptr[i]=bptr[r]; bptr[r]=tmptr;
    if (i>=k) r = i-1;
    if (i<=k) l = i+1;
  }
}

void sortbest (bptr, nbest, irelv)
struct beststr **bptr;
int nbest, irelv;
{
    int gap, i, j;
    struct beststr *tmp;

    for (gap = nbest/2; gap > 0; gap /= 2)
	for (i = gap; i < nbest; i++)
	    for (j = i - gap; j >= 0; j-= gap) {
	      if (bptr[j]->score[irelv] >= bptr[j + gap]->score[irelv]) break;
	      tmp = bptr[j];
	      bptr[j] = bptr[j + gap];
	      bptr[j + gap] = tmp;
	    }
}

void sortbestz(bptr, nbest)
struct beststr **bptr;
int nbest;
{
    int gap, i, j;
    struct beststr *tmp;

    for (gap = nbest/2; gap > 0; gap /= 2)
	for (i = gap; i < nbest; i++)
	    for (j = i - gap; j >= 0; j-= gap) {
	      if (bptr[j]->zscore >= bptr[j + gap]->zscore) break;
	      tmp = bptr[j];
	      bptr[j] = bptr[j + gap];
	      bptr[j + gap] = tmp;
	    }
}

void sortbeste(bptr, nbest)
struct beststr **bptr;
int nbest;
{
    int gap, i, j;
    struct beststr *tmp;

    for (gap = nbest/2; gap > 0; gap /= 2)
	for (i = gap; i < nbest; i++)
	    for (j = i - gap; j >= 0; j-= gap) {
	      if (bptr[j]->escore <= bptr[j + gap]->escore) break;
	      tmp = bptr[j];
	      bptr[j] = bptr[j + gap];
	      bptr[j + gap] = tmp;
	    }
}

extern double zs_to_Ec(float zs, long entries);

/*
extern float ks_dev;
extern int ks_df;
*/
extern char hstring1[];

void
prhist(FILE *fd, struct mngmsg m_msg,
       struct pstruct pst, 
       struct hist_str hist, 
       int nstats,
       struct db_str ntt,
       char *gstring2)
{
  int i,j,hl,hll, el, ell, ev;
  char hline[80], pch, *bp;
  int mh1, mht;
  int maxval, maxvalt, dotsiz, ddotsiz,doinset;
  float cur_e, prev_e, f_int;
  float max_dev, x_tmp;
  double db_tt;
  int n_chi_sq, cum_hl, max_i;


  fprintf(fd,"\n");
  
  if (pst.zsflag < 0 || nstats <= 10) {
    fprintf(fd, "%7ld residues in %5ld sequences\n", ntt.length,ntt.entries);
    fprintf(fd,"\n%s\n",gstring2);
    return;
  }

  max_dev = 0.0;
  mh1 = hist.maxh-1;
  mht = (3*hist.maxh-3)/4 - 1;

  if (!m_msg.nohist && mh1 > 0) {
    for (i=0,maxval=0,maxvalt=0; i<hist.maxh; i++) {
      if (hist.hist_a[i] > maxval) maxval = hist.hist_a[i];
      if (i >= mht &&  hist.hist_a[i]>maxvalt) maxvalt = hist.hist_a[i];
    }
    n_chi_sq = 0;
    cum_hl = -hist.hist_a[0];
    dotsiz = (maxval-1)/60+1;
    ddotsiz = (maxvalt-1)/50+1;
    doinset = (ddotsiz < dotsiz && dotsiz > 2);

    if (pst.zsflag>=0)
      fprintf(fd,"       opt      E()\n");
    else 
      fprintf(fd,"     opt\n");

    prev_e =  zs_to_Ec((float)(hist.min_hist-hist.histint/2),hist.entries);
    for (i=0; i<=mh1; i++) {
      pch = (i==mh1) ? '>' : ' ';
      pch = (i==0) ? '<' : pch;
      hll = hl = hist.hist_a[i];
      if (pst.zsflag>=0) {
	cum_hl += hl;
	f_int = (float)(i*hist.histint+hist.min_hist)+(float)hist.histint/2.0;
	cur_e = (float)zs_to_Ec(f_int,hist.entries);
	ev = el = ell = (int)(cur_e - prev_e + 0.5);
	if (hl > 0  && i > 5 && i < (90-hist.min_hist)/hist.histint) {
	  x_tmp  = fabs(cum_hl - cur_e);
	  if ( x_tmp > max_dev) {
	    max_dev = x_tmp;
	    max_i = i;
	  }
	  n_chi_sq++;
	}
	if ((el=(el+dotsiz-1)/dotsiz) > 60) el = 60;
	if ((ell=(ell+ddotsiz-1)/ddotsiz) > 40) ell = 40;
	fprintf(fd,"%c%3d %5d %5d:",
		pch,(i<mh1)?(i)*hist.histint+hist.min_hist :
		mh1*hist.histint+hist.min_hist,hl,ev);
      }
      else fprintf(fd,"%c%3d %5d :",
		   pch,(i<mh1)?(i)*hist.histint+hist.min_hist :
		   mh1*hist.histint+hist.min_hist,hl);

      if ((hl=(hl+dotsiz-1)/dotsiz) > 60) hl = 60;
      if ((hll=(hll+ddotsiz-1)/ddotsiz) > 40) hll = 40;
      for (j=0; j<hl; j++) hline[j]='='; 
      if (pst.zsflag>=0) {
	if (el <= hl ) {
	  if (el > 0) hline[el-1]='*';
	  hline[hl]='\0';
	}
	else {
	  for (j = hl; j < el; j++) hline[j]=' ';
	  hline[el-1]='*';
	  hline[hl=el]='\0';
	}
      }
      else hline[hl] = 0;
      if (i==1) {
	for (j=hl; j<10; j++) hline[j]=' ';
	sprintf(&hline[10]," one = represents %d library sequences",dotsiz);
      }
      if (doinset && i == mht-2) {
	for (j = hl; j < 10; j++) hline[j]=' ';
	sprintf(&hline[10]," inset = represents %d library sequences",ddotsiz);
      }
      if (i >= mht&& doinset ) {
	for (j = hl; j < 10; j++) hline[j]=' ';
	hline[10]=':';
	for (j = 11; j<11+hll; j++) hline[j]='=';
	hline[11+hll]='\0';
	if (pst.zsflag>=0) {
	  if (ell <= hll) hline[10+ell]='*';
	  else {
	    for (j = 11+hll; j < 10+ell; j++) hline[j]=' ';
	    hline[10+ell] = '*';
	    hline[11+ell] = '\0';
	  }
	}
      }

      fprintf(fd,"%s\n",hline);
      prev_e = cur_e;
    }
  }

  if (ntt.carry==0) {
    fprintf(fd, "%7ld residues in %5ld sequences\n", ntt.length, ntt.entries);
  }
  else {
    db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
    fprintf(fd, "%.0f residues in %5ld library sequences\n", db_tt, ntt.entries);
  }

  if (pst.zsflag>=0) {
    if (MAXSTATS < hist.entries)
      fprintf(fd," statistics extrapolated from %d to %ld sequences\n",
	      MAXSTATS,hist.entries);
    /*    summ_stats(stat_info); */
    fprintf(fd," %s\n",hist.stat_info);
    if (!m_msg.nohist && cum_hl > 0)
      fprintf(fd," Kolmogorov-Smirnov  statistic: %6.4f (N=%d) at %3d\n",
	      max_dev/(float)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
    if (m_msg.markx==10) {
      while ((bp=strchr(hist.stat_info,'\n'))!=NULL) *bp=' ';
      if (cum_hl <= 0) cum_hl = -1;
      sprintf(hstring1,"; mp_extrap: %d %ld\n; mp_stats: %s\n; mp_KS: %6.4f (N=%d) at %3d\n",
	      MAXSTATS,hist.entries,hist.stat_info,max_dev/(float)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
    }
  }
  fprintf(fd,"\n%s\n",gstring2);
  fflush(fd);
}

extern char prog_name[], *verstr;

void s_abort (char *p,  char *p1)
{
  int i;

  fprintf (stderr, "\n***[%s] %s%s***\n", prog_name, p, p1);
#ifdef PCOMPLIB
#ifdef PVM_SRC
  for (i=FIRSTNODE; i< nnodes; i++) pvm_kill(pinums[i]);
  pvm_exit();
#endif
#ifdef MPI_SRC
  MPI_Abort(MPI_COMM_WORLD,1);
  MPI_Finalize();
#endif
#endif
  exit (1);
}

#ifndef MPI_SRC
void w_abort (p, p1)
char *p, *p1;
{
  fprintf (stderr, "\n***[%s] %s%s***\n\n", prog_name, p, p1);
  exit (1);
}
#endif

#ifndef PCOMPLIB
/* copies from from to to shuffling */

extern int nrand(int);

void
shuffle(unsigned char *from, unsigned char *to, int n)
{
  int i,j; unsigned char tmp;

  if (from != to) memcpy((void *)to,(void *)from,n);

  for (i=n; i>1; i--) {
    j = nrand(i);
    tmp = to[j];
    to[j] = to[i-1];
    to[i-1] = tmp;
  }
  to[n] = 0;
}

/* copies from from to from shuffling, ieven changed for threads */
void
wshuffle(unsigned char *from, unsigned char *to, int n, int wsiz, int *ieven)
{
  int i,j, k, mm; 
  unsigned char tmp, *top;

  memcpy((void *)to,(void *)from,n);
	
  mm = n%wsiz;

  if (*ieven) {
    for (k=0; k<(n-wsiz); k += wsiz) {
      top = &to[k];
      for (i=wsiz; i>1; i--) {
	j = nrand(i);
	tmp = top[j];
	top[j] = top[i-1];
	top[i-1] = tmp;
      }
    }
    top = &to[n-mm];
    for (i=mm; i>1; i--) {
      j = nrand(i);
      tmp = top[j];
      top[j] = top[i-1];
      top[i-1] = tmp;
    }
    *ieven = 0;
  }
  else {
    for (k=n; k>=wsiz; k -= wsiz) {
      top = &to[k-wsiz];
      for (i=wsiz; i>1; i--) {
	j = nrand(i);
	tmp = top[j];
	top[j] = top[i-1];
	top[i-1] = tmp;
      }
    }
    top = &to[0];
    for (i=mm; i>1; i--) {
      j = nrand(i);
      tmp = top[j];
      top[j] = top[i-1];
      top[i-1] = tmp;
    }
    *ieven = 1;
  }
  to[n] = 0;
}

#endif

int
sfn_cmp(int *q, int *s)
{
  if (*q == *s) return *q;
  while (*q && *s) {
    if (*q == *s) return *q;
    else if (*q < *s) q++;
    else if (*q > *s) s++;
  }
  return 0;
}
